/*
 *  brick.h
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 01/04/09.
 *  Modified by Sergio Conde 16/09/2010
 */

#ifndef _brick_h_
#define _brick_h_

#include "component.h"
#include "commonutil.h"
#include "MATERIALS/material.h"

using namespace std;
// using namespace Eigen;

class Brick : public Component
{
public:
	Material *material; 		// Material data
// 	double mu,lambda;
	double detJ0;			
	double detF;
	
	MatrixXd gau;           // Gauss points matrix    (xi, eta, zeta) 8x3
	MatrixXd lob;			// Lobatto points matrix    (xi, eta, zeta) 8x3 or 27x3
	MatrixXd Q;         	// Coordinates (deformed) 8x3
	MatrixXd Qref;      	// Coordinates (reference) 8x3 
	Matrix3d J0, J0i;   	// Jacobian transformation between Euclidean coordinates and natural coordinates.
	MatrixXd H;             // Spatial derivatives shape function matrix
	MatrixXd N;             // Shape functions Matrix
	MatrixXd Nt;            // Shape functions Matrix for temperature
	MatrixXd Mt;
	VectorXd thetaVec;
	VectorXd thetaCVec;
	
	Matrix3d F;             // Deformation gradient
	Matrix3d CMat;			// Right Cauchy-Green Deformation Tensor: Matrix Format
	VectorXd C;             // Right Cauchy-Green Deformation Tensor: Vector Format
	MatrixXd vCConv;        // Vector of vectors, as many as Gauss points. Each matrix is 3x3
	VectorXd CMid;		// Right Cauchy-Green Deformation Tensor at mid-point
    VectorXd U;         // Displacement vector 24x1
	VectorXd the;     	// Vector of displacements derivatives 9x1
	MatrixXd vG;	    // Vector of matrices, as many as Gauss points. Each matrix is 8x3
	MatrixXd G;         // 9 x 24
	MatrixXd A;         // Matrix for Bnl calculation 6x9
	MatrixXd H0;        // Matrix for Bnl calculation 6x9
	MatrixXd B;         // Matrix Bnl (total) 6x24
	MatrixXd vBConv;    // Vector of matrices, as many as Gauss points. Each matrix is 8x3
	
	VectorXd vWConv;			// Vector of doubles, as many as Gauss points. Each matrix is 3x3
    double alpha;               // Constistent approximation factor for S
	double DeltaCDeltaC;		// DeltaC : DeltaC
	double SMidDeltaC;          // SMid : DeltaC
	double partial_qddotOld;
	double stEner;              //Strain Energy
	VectorXd S;       // Second Piola-Kirchhoff stress tensor (Vector format)
	MatrixXd vSConv;  // Vector of vectors, as many as Gauss points. Each matrix is 3x3
	MatrixXd vSMidConv;
	VectorXd SMid;    // Second Piola-Kirchhoff stress tensor at mid-ponit (Vector format)
	Matrix3d SMat;	  // Second Piola-Kirchhoff stress tensor (Matrix format)
	MatrixXd Sgor;	  // Stress dependent matrix for geometrical stiffness 9x9
	
    bool viscoMode;              //Viscoelastic internal variables integration type.
	int nInt;                   // Number of internal variables
	double beta;                // Constistent approximation factor for Qv
	double delta;               //
	double DeltaGamDeltaGam;	// DeltaGam : DeltaGam
	double QvMidDeltaGam;		// QvMid : DeltaGam
	VectorXd Gam;			// Internal Strain-like tensor. Vector of vector, as many as Gauss points.
	MatrixXd vGamConv;		// Converged Internal Strain-like tensor. Vector of vector, as many as Gauss points.
	VectorXd GamMid;		// Internal Strain-like tensor at mid-point. Vector of vector, as many as Gauss points.
	VectorXd auxGam;		//	Two internal variables
	VectorXd auxQv;			//	Two internal variables
	MatrixXd Kiv;		    // Tangent matrix for flow equation
	VectorXd resGam;		// Residue vector for flow equation
	double dissEner;        // Dissipated Energy
	VectorXd Qv;	  //Non-Equilibrium Stresses Tensor (Vector Format), as many as number of internal variables.
	MatrixXd vQvConv;  //Converged Non-Equilibrium Stresses Tensor (Vector Format), as many as number of internal variables.
	VectorXd QvMid;   //Non-Equilibrium Stresses Tensor at Midpoint (Vector Format), as many as number of internal variables.
	VectorXd Hv;	  //Hitory Term Tensor (Vector Format), as many as number of internal variables.
	VectorXd vTau;
	VectorXd vBeta;
	VectorXd vEta;
	VectorXd vMu;

	VectorXd integShape;

	MatrixXd I; 		//Identity Tensor in Voigt Format 6x6
	MatrixXd Cross;		// Matrix for Kq conservative calculation.
	VectorXd CrossKg;
	
	VectorXd Sigma;		// Cauchy stress tensor: Vector format
	Matrix3d SigmaMat;	// Cauchy stress tensor: Matrix format
	VectorXd PSigma;
	double VMStress;

	bool thermoMode;	// State variable: Temperature if thermoMode == true. Entropy if thermoMode == false
	double Theta;
	VectorXd vThetaConv;
	double Ent;
	double Entropy;
	VectorXd vEntConv;
	VectorXd vdEntConv;
	MatrixXd Kond;

	double kappa;
	VectorXd vkappaConv;

public:
	Brick();
	Brick(vector<string>&, string&, int);
	~Brick(){}

	/** Virtual functions from base */
	void setDofs(void);
	void initializeGlobalDofs(vector<double> &, vector<double> &, const double &);
	void calcFi(const vector<double>&);
	void calcFq(const int&);
	void calcFt(const double&, const double&, const int&);
	void calcKi(double&);
	void calcKq(int);
	void calcKdq(double&);
	void calcM(const vector<double>&, const double &);
	void calcMLumped(const vector<double>&, const double &);
	void setMatAtComp(Material*);
	string getMatName()	{return matName;}
	void initializeCompRes(const bool&, const double&, const int&, const bool);
//	void initializeCompTang(const bool&, const bool&);
	void initializeCompTang(const bool&);
	void calcEner(double&, double&, double&, double&, double&, double&, double&, int&, vector<double>&, vector<double>&);
	void calcEntropy(double&, double&);
	void calcAM(vector<double>&, vector<double>&, vector<double>&);
	void updateComp(void);
	void setInitialConditions(vector<double>&, vector<double>&){}
	double getTemperature(void){double T; return T;};
	void setIniDofs(vector<double>&);
	void calcPhi(const vector<double>&);
	void calcDPhi(const vector<double>&);
	void calcMSmooth(const vector<double>&, const double&);
	void calcFSmooth(const vector<double>&, const vector<double>&, const double&);
	int checkState(void);
	
	/** Local functions */
	void setQfromq();
	void makeU();	// U is the displacement vector
	void makeJ0(double &, double &, double &);
	void makevG();
	void makeG(int);
	void makeTheta();
	void makeA();
	void makeB();
	void makeF();
	void makeKq(int);
	void makeMLumped();
	void makeC();
	double calcShape(int&, double&, double&, double&);
	void calcIntegShape();
	double calcElementMass();
	
	void calcQv(int,int);	//Calculates the non-equilibrium stress in viscoelastic materials
	void calcHv(int,int);	//Calculates the history terms needed to calculte Qv.
    void calcQvConsistently(int); //Calculates the non-equilibrium stress as a partitioned discrete derivative of Psi respect to the Gamma.
	
	void integrateInternalStrains(int,int);
	void calcKivNumeric(int, int);
    //void calcKivConsistently(int, int);
	
	void integrateDamage(int,int);

	/** Friend functions */
};

#endif

